Analyse spatiale et territoriale

Formation Carthageo-Geoprisme 2023

Syllabus

Objectifs de la semaine

  • Développer les connaissances en analyse quantitative des données dans l’espace
  • Exploiter des données individuelles (demandes de valeurs foncières, recensement de la population)
  • Approfondir l’utilisation du logiciel R
  • Travailler en groupe et préparer une restitution commune
  • Répondre à une “commande” dans un temps limité

Organisation de la semaine

Planning des intervenants et contenu des journées

  • 09/10 (matin, Pierre Pistre) : Organisation de la semaine & présentation des données
  • 09/10 (après-midi, Claude Grasland) : Prise en main des données de recensement (fichier LOGEMT)
  • 10/10 (matin, Pierre Pistre) : Prise en main des données de valeurs foncières (DVF)
  • 10/10 (après-midi, Claude Grasland) : extractions des données pour la commande, et approfondissements statistiques & cartographiques
  • 11/10 : Accompagnement des projets de groupe (Hadrien Commenges)
  • 12/10 : Accompagnement des projets de groupe (Hadrien Commenges)
  • 13/10 : Finalisation des projets de groupe (matin) et restitution orale (après-midi)

Horaires et planning des salles (Campus des Grands Moulins, Université Paris Cité)

  • Horaires indicatifs : 9h30-12h30 et 13h30-17h.
  • Salles : du 09/10 au 11/10 : salle de cours 209 (Bâtiment Olympe de Gouges, 2ième étage) ; 12/10 : salle informatique 375 (Bâtiment Olympe de Gouges, 3ième étage) ; 13/10 : salle de cours 209 (Bâtiment Olympe de Gouges, 2ième étage)

Données utilisées

(sources) Base de données principales :

(sources) Fichiers complémentaires :

Exercice d’application en groupe

Périmètres de l’exercice

  • Objet d’étude : les dynamiques immobilières dans les principales agglomérations administratives de la “mégarégion parisienne” (hors Paris)
  • Population d’étude : les logements vendus récemment et leurs habitants
  • Espace(s) d’étude : les principales agglomérations ayant notamment le statut de “métrople” ou de “communauté urbaine” (hors Paris)
  • Echelles et mailles géographiques : département, commune, IRIS

Consignes de la “commande”

  • Contexte général : il est connu que Paris influence la dynamique des territoires bien au-delà de ses territoires les plus proches, à commencer par les principales agglomérations urbaines de la “mégarégion parisienne” (au sens élargi de : https://atlas-paris-mega-region.univ-rouen.fr/). Le secteur immobilier est particulièrement concerné du fait d’une progression sensible des prix des logements dans la plupart de ces territoires au cours à minima de la dernière décennie et d’un potentiel renforcement avec la crise COVID par une progression des installations résidentielles. Pour autant, la géographie des dynamiques immobilières internes à ces agglomérations urbaines reste mal connue de même que les facteurs explicatifs (influence parisienne, spécificités locales et régionales…) et les tendances récentes.
  • Commande : après avoir fait le constat de dynamiques immobilières souvent proches, les principales agglomérations urbaines de la “mégarégion parisienne” ont décidé de s’associer pour demander aux différentes antennes régionales de l’Insee la réalisation d’études rigoureuses et fines spatialement sur la situation immobilière de leurs différents territoires de compétence et périphéries proches. Le bilan de ces études sera présenté dans des publications synthétiques, intelligibles par le plus grand nombre et par les acteurs concernés, sous la forme de plusieurs 4 pages - comme produits régulièrement par l’Insee pour publier les résultats de ses études (par exemple, file:///home/pierre/T%C3%A9l%C3%A9chargements/ip1715.pdf ; file:///home/pierre/T%C3%A9l%C3%A9chargements/IR134_Notaires-IPLA_1T-23.pdf).

Modalités de l’exercice, restitution du travail et organisation durant la semaine

  1. Constitution de 6 (voire 7) groupes de 4-5 étudiant.e.s, mélangeant les profils de Master (= Lundi matin)
  2. Réflexion sur un angle problématique pour chaque cas d’étude, ainsi que les possibilités d’analyse et d’exploration des données : variables pertinentes, traitements envisagés… (= Lundi à Mercredi)
  3. Réalisation des analyses statistiques et cartographiques (= Mardi à Jeudi)
  4. Organisation, mise en forme des analyses et rédaction de la note écrite (= Mardi à Vendredi matin)
  5. Présentation orale (environ 10 minutes, sans autre support visuel que la note écrite) devant un “jury” composé des intervenants de la semaine (= Vendredi après-midi)

Cas d’étude

  • Caen (communauté urbaine)
  • Le Mans (communauté urbaine)
  • Orléans (métropole)
  • Reims (communauté urbaine)
  • Rouen (métropole)
  • Tours (métropole)

En option supplémentaire : Amiens (communauté d’agglomération)

Format du rendu (type 4 pages Insee)

  • Introduction et définition (thématique et statistique) de l’objet : thème et cas d’étude
  • Etat des lieux général des prix immobiliers dans le cas d’étude (données : DVF ; échelle : “agglomération” (rayon 40km) ; maille : commune)
  • Approfondissements sur les prix immobiliers par type de biens ou de localisation (données : DVF ; échelles : “agglomération” (rayon 40km) et zooms ; mailles : commune, IRIS)
  • Profils des acheteurs au sens des nouveaux arrivants (moins de 3 ans) qui sont propriétaires (données : recensement de la population ; échelles : “agglomération” (rayon 40km) et zooms ; mailles : commune, IRIS)

RP1 : Accquisition

Introduction

Ce programme permet de préparer l’ensemble des données statistiques et cartographiques utiles pour l’étude d’une agglomération urbaine à partir du fichier ménage du recensement de population 2019. Il charge de très gros fichiers contenus dans un dossier temporaire tmp et en extrait les informations utiles pour les placer dans un dossier data/nomagg . Il a besoin en entrée du code et du nom de la commune centre ainsi que du chemin vers le dossier de stockage des résultats

# Commune centre 
codectr<-  "35238"  # Code
namectr <- "Rennes"  # Nom

# Dossier de stockage
myrep <- "data/Rennes/"

# Choix du fichier ménage et de ses métadonnées
ficmen<-"tmp/FD_LOGEMTZC_2019.csv"
ficmen_meta <-"tmp/Varmod_LOGEMT_2019.csv"

Etape 1 : Données administratives

On commence par récupérer les informations adminstratives sur la zone d’étude

Chargement du fichier administatif

adm<-read.table("tmp/maps/table-appartenance-geo-communes-23.csv",
                sep=";", 
                header=T, fileEncoding = "UTF-8",
                quote = '"')
head(adm)

Identification des codes de l’aire urbaine et de l’EPCI

On extrait la ligne correspondant à notre commune centre et on en déduit le code de son epci et de son aire urbaine

adm_ctr<-adm %>% filter(adm$CODGEO==codectr)
adm_ctr

my_EPCI <- adm_ctr$EPCI
my_AA <- adm_ctr$AAV2020

Extraction des informations

On extrait toutes les donnes relatives à notre aire urbaine et on les sauvegarde.

AA_adm <- adm %>% filter(AAV2020 == my_AA)
saveRDS(AA_adm, paste0(myrep,"AA_adm.RDS"))

EPCI_adm <- adm %>% filter(EPCI == my_EPCI)
saveRDS(EPCI_adm, paste0(myrep,"EPCI_adm.RDS"))

CTR_adm <- adm %>% filter(CODGEO == codectr)
saveRDS(CTR_adm, paste0(myrep,"CTR_adm.RDS"))

Etape 2 : Données géométriques

On va maintenant extraire les données géométriques relatives au contour des communes et des IRIS

Communes

On charge le fichier des communes de France et on extrait uniquement celles de la zone d’étude.

## Communes de France
mapcom <- st_read("tmp/maps/COMMUNE.shp", quiet=T)

## Aire urbaine
AA_adm<-readRDS(paste0(myrep,"AA_adm.RDS"))
AA_mapcom <- mapcom %>% filter(INSEE_COM %in% AA_adm$CODGEO)
saveRDS(AA_mapcom,paste0(myrep,"AA_mapcom.RDS"))

## EPCI
EPCI_adm<-readRDS(paste0(myrep,"EPCI_adm.RDS"))
EPCI_mapcom <- mapcom %>% filter(INSEE_COM %in% EPCI_adm$CODGEO)
saveRDS(EPCI_mapcom,paste0(myrep,"EPCI_mapcom.RDS"))

## Commune centre
CTR_adm<-readRDS(paste0(myrep,"CTR_adm.RDS"))
CTR_mapcom <- mapcom %>% filter(INSEE_COM %in% CTR_adm$CODGEO)
saveRDS(CTR_mapcom,paste0(myrep,"CTR_mapcom.RDS"))
AA_mapcom<-readRDS(paste0(myrep,"AA_mapcom.RDS"))
EPCI_mapcom<-readRDS(paste0(myrep,"EPCI_mapcom.RDS"))
CTR_mapcom<-readRDS(paste0(myrep,"CTR_mapcom.RDS"))
plot(AA_mapcom$geometry, col="lightyellow", main="zone d'étude")
plot(EPCI_mapcom$geometry,col="orange", add=T)
plot(CTR_mapcom$geometry,col="red", add=T)

IRIS

On ne procède à l’extraction des IRIS que pour l’EPCI et la commune centre

## Communes de France
mapiris <- st_read("tmp/maps/CONTOURS-IRIS.shp", quiet=T)


## EPCI
EPCI_adm<-readRDS(paste0(myrep,"EPCI_adm.RDS"))
EPCI_mapiris <- mapiris %>% filter(INSEE_COM %in% EPCI_adm$CODGEO)
saveRDS(EPCI_mapiris,paste0(myrep,"EPCI_mapiris.RDS"))

## Commune centre
CTR_adm<-readRDS(paste0(myrep,"CTR_adm.RDS"))
CTR_mapiris <- mapiris %>% filter(INSEE_COM %in% CTR_adm$CODGEO)
saveRDS(CTR_mapiris,paste0(myrep,"CTR_mapiris.RDS"))
EPCI_mapiris<-readRDS(paste0(myrep,"EPCI_mapiris.RDS"))
CTR_mapiris<-readRDS(paste0(myrep,"CTR_mapiris.RDS"))
plot(EPCI_mapiris$geometry, col="lightyellow", main="EPCI",border="gray50", lwd=0.4)
plot(CTR_mapiris$geometry,col="orange", add=T, border="gray50", lwd=0.4)
plot(EPCI_mapcom$geometry,col=NA, border="black",lwd=1,add=T)

Etape 3 : logements ordinaires en 2019

Nous partirons des fichiers détail de l’INSEE car, à la différence des tableaux prédéfinis, ils permettent virtuellement toutes les formes de croisement d’indicateurs. Ils sont évidemment beaucoup plus volumineux, mais ce sera justement l’occasion pour les étudiants en data mining d’être confrontés à des problèmes d’optimisation et de big data. On trouve leur description détaillée sur le site de l’INSEE

Compte-tenu de la taille des fichiers, nous travaillerons sur une partition de la France en zones

  • Zone A : Région Île-de-France (région 11) ;
  • Zone B : Régions Centre-Val de Loire (région 24), Bourgogne-Franche-Comté (région 27), Normandie (région 28) et Hauts-de-France (région 32) ;
  • Zone C : Régions Grand Est (région 44), Pays de la Loire (région 52) et Bretagne (région 53);
  • Zone D : Régions Nouvelle-Aquitaine (région 75) et Occitanie (région 76) ;
  • Zone E : Régions Auvergne-Rhônes-Alpes (région 84), Provence-Alpes-Côte d’Azur (région 93), Corse (région 94), Guadeloupe (région 01), Martinique (région 02), Guyane (région 03) et La Réunion (région 04).

Lecture du fichier des ménages

Nous commençons par lire le fichiers de la zone qui nous intéresse (au format .csv). Nous utilisons pour cela la fonction fread du package data.table qui est très rapide mais qui commet une erreur sur le code communal qu’on doit corriger

# Lit le fichier
men<-fread(ficmen)


# Ajoute la variable INSEE_COM en format caractère en rétablissant le zéro
code <-as.character(men$COMMUNE)
code[nchar(code)==4]<-paste0("0",code[nchar(code)==4])
men$INSEE_COM<-code

Extraction des ménages de la zone d’étude

## Aire urbaine
AA_adm<-readRDS(paste0(myrep,"AA_adm.RDS"))
AA_men <- men %>% filter(INSEE_COM %in% AA_adm$CODGEO)
saveRDS(AA_men,paste0(myrep,"AA_men.RDS"))

## EPCI
EPCI_adm<-readRDS(paste0(myrep,"EPCI_adm.RDS"))
EPCI_men <- men %>% filter(INSEE_COM %in% EPCI_adm$CODGEO)
saveRDS(EPCI_men,paste0(myrep,"EPCI_men.RDS"))

## Commune centre
CTR_adm<-readRDS(paste0(myrep,"CTR_adm.RDS"))
CTR_men <- men %>% filter(INSEE_COM %in% CTR_adm$CODGEO)
saveRDS(CTR_men,paste0(myrep,"CTR_men.RDS"))

Chargement des méta-données

Il ne faut surtout pas oublier le fichier des métadonnées qui va permettre de recoder facilement tous les facteurs et de décoder les chiffres correspondant aux classes. On va donc le transformer au format R puis l’enregistrer également dans le dossier data. On élimine toutefois les codes des variables COMMUNES, IRIS et TRIRIS qui prennent trop de place.

# Lecture du fichier de métadonnées
meta<-fread(ficmen_meta)

# Elimine les modalités IRIS, COMMUNE et TRIRIS
meta <-meta %>% filter(!(COD_VAR %in% c("COMMUNE", "IRIS","TRIRIS")))

# Enregistrement dans le dossier data
saveRDS(object = meta,
        file = paste0(myrep,"men_meta.RDS"))

RP2 : Agrégation

Le format sf (spatial features)

La cartographie et plus généralement les opérations géométriques sur des données spatiales dans R peuvent facilement être effectuées avec le package sf (spatial features) qui crée des objets uniques rassemblant à la fois

  • un tableau de données (l’équivalent du fichier .dbf)
  • une géométrie (l’équivalent du fichier .shp)
  • une projection (l’équivalent du fichier .prj)

Lorsqu’on récupère des fonds de carte au format shapefile (.shp) ou dans d’autres formats standards comme GeoJson, la première tâche consiste donc à les convertir au formt sf afin de pouvoir les utiliser facilement dans R. L’importation se fait à l’aide de l’instruction st_read en indiquant juste le nom du fichier .shp à charger. Les autres fichiers (.dbf ou .proj) seront lus également et intégrés dans l’objet qui hérite de la double classe data.frame et sf.

Etapes de préparation des données

Dans notre exemple, nous allons suivre les étapes suivantes :

  1. Préparer les données statistiques par IRIS dans un data.frame
  2. Charger un fonds de carte par IRIS au format sf
  3. Effectuer une jointure entre les deux fichiers par le code IRIS
  4. Sauvegarder le résultat
  5. Agréger les données statistiques et géométriques par commune
  6. Sauvegarder le résultat.

Préparer les données statistiques

On importe le fichier des individus et on corrige le code IRIS des communes non divisées en IRIS.

myrep<-"data/Rennes/"
tab_ind<-readRDS(paste0(myrep,"EPCI_men.RDS"))
tab_ind$IRIS[tab_ind$IRIS=="ZZZZZZZZZ"]<-
  paste0(tab_ind$INSEE_COM[tab_ind$IRIS=="ZZZZZZZZZ"],"0000")

head(tab_ind[,1:5],3)
FALSE    COMMUNE   ARM      IRIS ACHL AEMM
FALSE 1:   35001 ZZZZZ 350010101  B12 1976
FALSE 2:   35001 ZZZZZ 350010101  B12 1977
FALSE 3:   35001 ZZZZZ 350010101  B12    0

Selectionner les lignes et colonnes

On décide de ne sélectionner que les ménages de propriétaires installés depuis moins de 3 ans et on retient comme variable le type de logement que l’on recode en deux catégories

tab_ind2 <- tab_ind %>% filter(as.numeric(ANEM) < 3,        # durée d'occupation < 3 ans
                               STOCD =="10",                # propriétaire de son logement
                               TYPL %in% c(1,2)) %>%        #  Maison ou Appartement  
                        mutate(TYPL = as.factor(TYPL)) %>%
                        select(IPONDL, IRIS, TYPL) 
levels(tab_ind2$TYPL)<-c("Maison","Appt")

knitr::kable(head(tab_ind2,5),digits=2)
IPONDL IRIS TYPL
1.07 350010101 Maison
1.03 350010101 Maison
1.07 350010101 Maison
0.99 350010101 Maison
1.11 350010101 Maison

Agréger les données

On commence par créer un tableau long croisant les deux variables et leur effectif pondéré :

tab_long<- tab_ind2 %>%
           group_by(IRIS,TYPL)%>%
           summarise(NB=sum(IPONDL))


knitr::kable(head(tab_long,5),digits=2)
IRIS TYPL NB
350010101 Maison 36.13
350010101 Appt 9.71
350010102 Maison 106.53
350010102 Appt 57.60
350220000 Maison 21.89

Pivoter le tableau

Puis on fait “pivoter” le tableau pour l’obtenir en format large :

tab_large <- tab_long %>% pivot_wider(id_cols = IRIS, 
                                      names_from = TYPL,
                                      names_prefix = "T_",
                                      values_from = NB,
                                      values_fill = 0)

knitr::kable(head(tab_large,5),digits=2)
IRIS T_Maison T_Appt
350010101 36.13 9.71
350010102 106.53 57.60
350220000 21.89 0.99
350240101 142.85 41.34
350240102 34.58 10.50

Ajouter de nouvelles variables

On ajoute de nouvelles variables telles que le nombre total de nouveaux propriétaires et la part des propriétaires de maisons parmi eux.

tab<- tab_large %>% mutate(T_TOT = T_Maison+T_Appt,
                           Maison_pct = 100*T_Maison/T_TOT)

knitr::kable(head(tab,5),digits=c(0,0,0,0,1))
IRIS T_Maison T_Appt T_TOT Maison_pct
350010101 36 10 46 78.8
350010102 107 58 164 64.9
350220000 22 1 23 95.7
350240101 143 41 184 77.6
350240102 35 10 45 76.7

Examiner la distribution statistique

On examine l’histogramme donnant distribution statistique du % de nouveux propriétaires en maison individuelle.

programme

p <- ggplot(tab) + aes (x = Maison_pct) +
                   geom_histogram(breaks = c(0,10,20,30,40,50,
                                             60,70,80,90, 100)) +
                   scale_x_continuous("% de Maisons") +
                   scale_y_continuous("Nombre d'IRIS") +
                   ggtitle(label = "Type de logement des nouveaux propriétaires",
                           subtitle = "Source : INSEE, RP 2019 - Métropole de Rennes")
                            
p

Charger les données géométriques

On importe le fichier des iris de l’agglomération qui est de type sf (spatial features)

map_iris <- readRDS(paste0(myrep,"EPCI_mapiris.RDS"))
map_iris<-map_iris[,c(5,6,3,2)]
names(map_iris)<-c("IRIS","NOM_IRIS","COM","NOM_COM","geometry")

class(map_iris)
FALSE [1] "sf"         "data.frame"
knitr::kable(head(as.data.frame(map_iris)[,1:5],2))
IRIS NOM_IRIS COM NOM_COM geometry
352381204 Les Cloteaux Rennes 35238 MULTIPOLYGON (((350568 6786…
352750000 Saint-Gilles Saint-Gilles 35275 MULTIPOLYGON (((337255.9 67…

Visualisation du fonds iris avec sf

On peut facilement faire un plot de la colonne geometry du fichier sf

plot(map_iris$geometry,col="lightyellow")

Jointure des données IRIS et du fonds de carte

programme

map_iris_tab<-merge(map_iris,tab,
                   by.x="IRIS",by.y="IRIS",
                   all.x=T,all.y=F)

knitr::kable(head(map_iris_tab,3),digits=2)
IRIS NOM_IRIS COM NOM_COM T_Maison T_Appt T_TOT Maison_pct geometry
350010101 Ouest Acigné 35001 36.13 9.71 45.84 78.82 MULTIPOLYGON (((360989.4 67…
350010102 Est Acigné 35001 106.53 57.60 164.13 64.90 MULTIPOLYGON (((364373.3 67…
350220000 Bécherel Bécherel 35022 21.89 0.99 22.88 95.65 MULTIPOLYGON (((333530.4 68…

Cartographie rapide avec sf

Une astuce pour visualiser rapidement une carte choropléthe avec le seul packahe sf :

plot(map_iris_tab[,"Maison_pct"])

Sauvegarde du fichier par IRIS

Au format RDS (pour réutilisation dans R)

saveRDS(map_iris_tab,paste0(myrep,"map_iris_typl.RDS"))

Au format shapefile (pour Magrit ou Qgis)

st_write(map_iris_tab,paste0(myrep,"map_iris_typl.shp"))

Agrégation statistique + géométriques

Grâce aux nouveaux packages de R (dplyr et sf) il est possible d’agréger simultanément les statistiques et les géométries après les avoir stockés dans un même objet de type “sf”

Du coup, on peut gagner beaucoup de temps dans les traitements et les analyses cartographiques, en particulier si l’on veut tester différents niveaux d’agrégation.

Agrégation des IRIS en communes

L’agrégation est très facile et elle concerne à la fois les variables (de stock) et les geometries

map_com_tab <- map_iris_tab %>% 
  group_by(COM, NOM_COM) %>% 
  summarise(T_Maison=sum(T_Maison,na.rm=T), 
            T_Appt=sum(T_Appt,na.rm=T)) %>%
  st_cast("MULTIPOLYGON")

map_com_tab <- map_com_tab %>%  mutate(T_TOT = T_Maison+T_Appt,
                                  Maison_pct = 100*T_Maison/T_TOT) 

Agrégation des iris en communes

résultat statistique

knitr::kable(head(st_drop_geometry(map_com_tab),3),digits=c(0,0,0,0,0,1))
COM NOM_COM T_Maison T_Appt T_TOT Maison_pct
Acigné 35001 143 67 210 67.9
Betton 35024 291 82 373 78.0
Bourgbarré 35032 162 24 186 87.2

résultat géométrique

plot(map_com_tab$geometry,col ="lightyellow")

Examiner la distribution statistique

On examine l’histogramme donnant distribution statistique du % de nouveaux propriétaires en maison par commune.

p <- ggplot(map_com_tab) +aes (x = Maison_pct) +
                   geom_histogram(breaks = c(0,10,20,30,40,50,
                                             60,70,80,90, 100)) +
                   scale_x_continuous("% de Maisons") +
                   scale_y_continuous("Nombre de communes") +
                   ggtitle(label = "Type de logement des nouveaux propriétaires",
                           subtitle = "Source : INSEE, RP 2019, Agglomération de Rennes")
p                            

Cartographie rapide

plot(map_com_tab[,"Maison_pct"])


RP3 : Cartographie

Le package map_sf

Le package mapsf permet de réaliser des cartes statiques de très haute qualité. Il a en effet été mis au point par des cartographes et des géomaticiens professionnels de l’UMS RIATE. Il prend la suite du package cartography dont la maintenance demeurera assuré quelque temps encore mais ne fera plus l’objet de développements futurs. Le package mapsf présente l’avantage d’être totalement compatibvle avec le package sf ce qui n’était pas autant le cas pour le package cartography, plus ancien, et créé pour être compatible avec l’ancien package sp.

On trouvera la documentation du package mapsf à l’adresse suivante :

https://riatelab.github.io/mapsf/index.html

Création d’un template cartographique

Nous allons dans un premier temps apprendre à créer un fonds de carte vierge mais comportant tout l’habillage nécessaire (“template”). Pour cela nous allons charger différentes couches cartographiques correspondant respectivement au département, aux communes et aux iris :

myrep<-"data/Rennes/"
map_iris<-readRDS(paste0(myrep,"map_iris_typl.RDS"))
map_com<-readRDS(paste0(myrep,"map_com_typl.RDS"))

tracé d’un fonds de carte vierge

La fonction mf_map() avec le paramètre type = "base"permet de tracer une carte vide

 mf_map(map_iris, type = "base")

Superposition de couches

On peut toutefois ajouter toute une série de paramètres supplémentaire (col=, border=, lwd=, …) et superposer plusieurs fonds de carte avec le paramètre add = TRUE. L’ajout de la fonction layout permet de rajouter un cadre une légende.

# Trace les Iris avec des paramètres
mf_map(map_iris,  
       type = "base", 
       col = "lightyellow", 
       border="gray50",
       lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com,  
       type = "base", 
       col = NA,
       border="red",
       lwd=0.6,
       add = TRUE)

# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et IRIS de la zone d'étude", 
          credits = "Sources : IGN et INSEE")

Ajout d’un thème

On peut finalement modifier l’ensemble de la carte en lui ajoutant une instruction mf_theme() qui peut reprendre des styles existants ( “default”, “brutal”, “ink”, “dark”, “agolalight”, “candy”, “darkula”, “iceberg”, “green”, “nevermind”, “jsk”, “barcelona”) mais aussi créer ses propres thèmes

#Choix du thème
mf_theme("darkula")
# Trace les Iris avec des paramètres
mf_map(map_iris,  
       type = "base", 
       col = "lightyellow", 
       border="gray50",
       lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com,  
       type = "base", 
       col = NA,
       border="red",
       lwd=0.6,
       add = TRUE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et IRIS de la zone d'étude", 
          credits = "Sources : IGN et INSEE")

Ajout de texte

On peut ajouter une couche de texte avec la fonction mf_label(). Par exemple, on va ajouter à la carte précédente le nom des communes

programme

mf_theme("agolalight")

# Trace les Iris avec des paramètres
mf_map(map_iris,  
       type = "base", 
       col = "lightyellow", 
       border="gray50",
       lwd=0.3)

# Ajoute les contours des communes
mf_map(map_com,  
       type = "base", 
       col = NA,
       border="red",
       lwd=0.6,
       add = TRUE)

# Ajoute une couche de labels
mf_label(map_com, 
         var="COM",
         cex=0.4,
         col="blue",
         overlap = FALSE)

# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et IRIS de la zone d'étude", 
          credits = "Sources : IGN et INSEE")

Carte de stock

Une carte de stock représente la localisation de quantités que l’on peut aditionner et dont le total a un sens. Par exemple un nombre d’habitants, un nombre de ménages, un nombre d’automobiles. Ce quantités doivent être représentées par des figures (cercles, carrés, …) dont la surface est proportionelle au stock afin que l’oeil du lecteur puisse les aditionner visuellement.

Dans le package mapsf, on réalise ce type de carte à l’aide de la fonction mf_map()en lui donnant le paramètre type="prop".

On va tenter à titre d’exemple de représenter la distribution du nombre de nouveaux propriétaires.

Carte de stock minimale

Les instructions minimales sont les suivantes :

# Trace les contours des communes
mf_map(x= map_iris, 
       type = "base")

# Ajoute le nombre de ménages par IRIS
mf_map(x =map_iris, 
      type ="prop",
      var = "T_TOT",
      add=TRUE)

Mais le résultat est peu satisfaisant car les cercles sont trop grands. Il faut en pratique toujours effectuer un réglage de ceux-ci avec l’instruction inches=

Carte de stock habillée

mf_theme("agolalight")
mf_map(map_iris, type = "base",  
       col = "lightyellow",border="gray80", lwd=0.3)
mf_map(map_com, type = "base", 
       col = NA,border="black",lwd=1,add = TRUE)

mf_map(map_iris, var = "T_TOT",type = "prop",
  inches = 0.05, col = "red",leg_pos = "left",  
  leg_title = "Nombre de ménages", add=TRUE)

mf_layout(title = "Distribution des nouveaux propriétaires", 
          frame = TRUE,
          credits = "Sources : IGN et INSEE")

Carte choroplèthe

Une carte choroplèthe ou d’intensité représente un phénomène relatif dont la somme n’a pas de sens. Par exemple, il serait absurde d’aditionner les % de logement HLM des IRIS du Val de Marne. Ces variables d’intensité caractèrisent donc l’état général d’une zone (choros) et elles vont être représentées par une couleur appliquée à toute la surface de la zone, d’où leur nom de cartes choroplèthes.

La fonction du package mapsf adaptée aux variables d’intensité est la fonction mf_map()munie du paramètre type = "choro".

On va prendre l’exemple du % de maisons parmi les nouveaux propriétaires

Carte choroplèthe minimale

Si on ne précise rien, la carte est réalisée à l’aide de la palette par défaut avec un découpage des classes en quantiles (effectifs égaux).

# Carte choroplèthe
mf_map(
  x = map_iris, 
  var = "Maison_pct",
  type = "choro")

Carte choroplèthe habillée

On peut arriver à une carte beaucoup plus satisfaisante en contrôlant l’ensemble des paramètres de couleur et de découpage des classes. Puis en superposant les contours de communes au dessus de la carte des IRIS pour faciliter le repérage.

mybreaks = c(0, 10,20,30,40,50,
             60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5), 
                    pal = c("Greens", "Reds"))
# Carte choroplèthe des iris
mf_map( map_iris, var = "Maison_pct",
        type = "choro",
        breaks = mybreaks,pal = mypal, 
        border=NA,
       col_na = "gray80",leg_title = "% maisons",
       leg_val_rnd = 0)
# Contour des communes et cadre
mf_map(map_com, type = "base", col = NA,
       border="black",lwd=1,add = TRUE)
mf_layout(title = "Les nouveaux propriétaires de maisons", 
          frame = TRUE,
          credits = "Sources : IGN et INSEE")

Carte stock + choroplèthe (1)

On peut combiner les deux modes cartographiques par superposition :

mf_theme("agolalight")

# Choisit les classes
mybreaks = c(0,20,40,60,80,100)

# Trace la carte choroplèthe
mf_map(
  x = map_iris, 
  var = "Maison_pct",
  breaks = mybreaks,
 # pal=mypal,
 type = "choro",
  border="white",
  col_na = "gray80",
 lwd=0.3,
 leg_title = "% Maisons", 
 leg_val_rnd = 0,
  
)

# Ajoute les cercles proportionnels

mf_map(
  x =map_iris, 
  var = "T_Maison",
  type = "prop",
  inches = 0.06, 
  col = "red",
  leg_pos = "right",  
  leg_title = "Effectif",
  add=TRUE
)
# Ajoute les contours des communes
mf_map(map_com, 
       type = "base", 
       col = NA,
       border="black",
       lwd=1,
       add = TRUE)

# Ajoute un cadre, un titre et des sources
mf_layout(title = "Les nouveaux propriétaires de maison", 
          frame = TRUE,
          credits = "Sources : IGN et INSEE")

Carte stock + choroplèthe (2)

Mais les cercles dissimuent alors les plages de couleur, aussi on peut utiliser le type prop_choro qui place la variable choroplèthe à l’intérieur des cercles

mf_theme("agolalight")
mybreaks = c(0, 10,20,30,40,50,60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5), pal = c("Greens","Reds"))
mf_map(map_iris, type = "base",  
       col = "gray80",border="white", lwd=0.3)
mf_map(map_com, type = "base", 
       col = NA,border="white",lwd=1,add = TRUE)
mf_prop_choro( x = map_iris,  var = c("T_TOT", "Maison_pct"), 
  inches = 0.06, col_na = "grey", pal=mypal,
  breaks = mybreaks, nbreaks = 4, lwd = 0.1,
  leg_pos = c("right", "left"),leg_val_rnd = c(0,0),
  leg_title = c("Nouveaux propriétaires", "% Maison"),
  add = TRUE)
mf_layout(title = "Les nouveaux propriétaires de maisons", 
          frame = TRUE, credits = "Sources : IGN et INSEE")


DVF1 : Acquisition

Introduction

Objectif

On se propose de collecter l’ensemble des ventes de maison ou d’appartement dans un rayon de 50 km autour d’une ville. On choisit ici Rennes en prenant la place du Parlement de Bretagne comme centre de référence

Paramètres

# Commune centre 
codectr<-  "35238"  # Code
namectr <- "Rennes"  # Nom
latctr <-  -1.67789
lonctr <-  48.11204
# Choix du rayon de collecte des dvf (en mètres)
rayon <- 50000
# Dossier de stockage
myrep <- "data/Rennes/"

Données dvf

Les données sur les dvf sont disponibles sur le site publicopendatasof:

https://public.opendatasoft.com

On s’intéresse plus particulièrement aux données dvf géolocalisées accessibles ici :

https://public.opendatasoft.com/explore/dataset/buildingref-france-demande-de-valeurs-foncieres-geolocalisee-millesime/

Sélection des informations

Le fichier comporte 28 millions d’enregistrement ce qui est évidemment beaucoup … On va donc procéder à une sélection en se servant des différents onglets disponibles.

A titre d’exemple, nous allons essayer de télécharger les enregistrements vérifiant les conditions suivantes :

  • localisation dans un rayon de 60 km autour d’Amiens
  • type de transaction : Vente
  • type de bien : Maison ou appartement

Nous constatons qu’il y a 133294 enregistrements vérifiant ces conditions.

Lien de téléchargement

Nous pouvons récupérer les données soit au format .csv pour les analyse statistiques, soit au format .geojson pour les analyses spatiales. Mais dans ce cas on ne peut pas dépasser la taille de 50 000 enregistrement. Il vaut donc mieux télécharger au format .csv puisque ce fichier comporte les latitudes et longitudes ce qui nous permettra de créer un fichier cartographique a posteriori

On peut effectuer le téléchargement depuis le site web ou bien juste enregistrer le lien de téléchargement et effectuer l’opération dans R.

Pour trouver le lien on passe la souris au dessus du lien “télécharger les données au format .csv / seulement les enregistrements sélectionnées” et on effectue un click droit suivi de “copier le lien

Construction d’une API

On récupère le lien de téléchargement du fichier .csv :

https://public.opendatasoft.com/api/explore/v2.1/catalog/datasets/buildingref-france-demande-de-valeurs-foncieres-geolocalisee-millesime/exports/csv?lang=fr&refine=nature_mutation%3A%22Vente%22&refine=type_local%3A%22Maison%22&refine=type_local%3A%22Appartement%22&facet=facet(name%3D%22nature_mutation%22%2C%20disjunctive%3Dtrue)&facet=facet(name%3D%22type_local%22%2C%20disjunctive%3Dtrue)&where=(distance(%60geo_point%60%2C%20geom%27POINT(2.24771118722856%2049.91161703449722)%27%2C%2060298.74879666418m))&timezone=Europe%2FParis&use_labels=true&delimiter=%3B

A première vue ceci paraît assez obscur mais on réalise assez vite que cette URL reprend l’ensemble des spécifications envoyées pour sélectionner nos enregistrements. On peut souligner en gras les parties utiles :

  • &refine=nature_mutation%3A%22Vente%22&
  • &refine=type_local%3A%22Maison%22
  • &refine=type_local%3A%22Appartement%22
  • &where=(distance(%60geo_point%60%2C%20geom%27
  • POINT(2.26968%2049.90807)%27%2C%2060037m%20))

Nous allons donc créer un lien qui extrait automatiquement les dvf en modifiant juste la position du point central et le rayon. On prend cette fois-ci nos paramètres de Rennes

myurl<-paste0("https://public.opendatasoft.com/api/explore/v2.1/catalog/datasets/buildingref-france-demande-de-valeurs-foncieres-geolocalisee-millesime/exports/csv?lang=fr&refine=nature_mutation%3A%22Vente%22&refine=type_local%3A%22Maison%22&refine=type_local%3A%22Appartement%22&facet=facet(name%3D%22nature_mutation%22%2C%20disjunctive%3Dtrue)&facet=facet(name%3D%22type_local%22%2C%20disjunctive%3Dtrue)&where=(distance(%60geo_point%60%2C%20geom%27POINT(", latctr,"%20",lonctr,")%27%2C%20",rayon,"m))&timezone=Europe%2FParis&use_labels=true&delimiter=%3B")

Téléchargement

On lance la requête avec l’instruction download.file() et on stocke le résultat dans notre répertoire de travail :

download.file(url = myurl, 
              destfile = paste0(myrep,"dvf.csv"))

essai de l’URL ‘https://public.opendatasoft.com/…’ downloaded 38.4 MB

L’opération s’execute en un peu moins d’une minute dans cet exemple. Mais il peut arriver qu’elle échoue lorsque le serveur est saturé et dans ce cas il faut recommencer.

DVF2 : Nettoyage

Introduction

Objectif

On se propose de nettoyer les donénes relative à l’ensemble des ventes de maison ou d’appartement dans un rayon de 50 km autour d’une ville. On rappelle les paramètres qui ont été utilisés pour acquérir ces données.

Paramètres

# Commune centre 
codectr<-  "35238"  # Code
namectr <- "Rennes"  # Nom
latctr <-  -1.67789
lonctr <-  48.11204
# Choix du rayon de collecte des dvf (en mètres)
rayon <- 50000
# Dossier de stockage
myrep <- "data/Rennes/"
# Fichier dvf brut
myfic <-"dvf.csv"

Merci à Boris et Florent !

Nous allons effectuer le nettoyage en nous appuyant sur les propositions de Boris Mericskay et Florent Demoraes parues dans https://journals.openedition.org/cybergeo/39583

Les auteurs fournissent en effet un accès à tous les programmes utilisés dans leur article via une archive github : https://github.com/ESO-Rennes/Analyse-Donnees-DVF

Chargement du jeu de données

On adapte le programme de Boris et Florent en utilisant la fonction fread() du package data.table car elle est plus rapide/

DVF<-fread(paste0(myrep,myfic))

Etape 1 : Sélection des mutations de type “Ventes” de “Maison” et “Appartement’

Programme

etape1 <- DVF %>% filter(`Nature de la mutation` == "Vente")
etape1bis <- etape1 %>% filter(`Type de local` == "Maison" |
                              `Type de local` == "Appartement")

Remarque

Etape inutile si on a déjà fait ce double filtrage à l’aide de l’API

Etape 2 : Sélection et renommage des variables

etape2 <- etape1bis %>% 
          select(id = `Identifiant de mutation (Etalab)`,
                 disposition = `Numéro de disposition`  , 
                 parcelle= `Identifiant de la parcelle cadastrale`,
                 date = `Date de la mutation`,
                 nature = `Nature de la mutation`, 
                 INSEE_COM = `Code INSEE de la commune`,
                 INSEE_DEP =  `Code INSEE du département` , 
                 type = `Type de local` ,
                 surface = `Surface réelle du bâti` , 
                 piece = `Nombre de pièces principales`, 
                 prix = `Valeur foncière`  , 
                 latitude = `Latitude`, 
                 longitude = `Longitude`)

Etape 3 : Remplacement des cellules vides par des NA et suppression des NA

etape2[etape2==""] <- NA
etape3 <- etape2 %>% na.omit()

Etape 4 : Sélections des mutations simples

Regrouper les transactions selon l’ID, la surface et la vente

unique <- etape3 %>% distinct(id, prix, surface)
nbunique <- unique %>% group_by(id) %>% summarise(nb = n())

Suppression des doublons et des mutations multiventes

etape4 <- nbunique %>% filter(nb==1)

Etape 5 : Jointure attributaire pour récupérer les informations de la mutation

jointure

merge <- merge(etape4,etape3, by="id")
etape5 <- merge %>% distinct(id, .keep_all = TRUE) %>% 
select(id, date, type, nature, INSEE_COM, prix,
      surface, piece, latitude, longitude)

Modification des formats des colonnes

etape5$prix <- as.numeric(etape5$prix)
etape5$surface <- as.numeric(etape5$surface)
etape5$piece <- as.numeric(etape5$piece)

Suppression des valeurs aberrantes

Fixer un seuil minimal et maximal des prix (percentiles)

quantile(etape5$prix, c(0.01, 0.99))
    1%    99% 
 17000 616320 

Fixer un seuil minimal et maximal des prix (histogramme)

hist(etape5$prix, breaks = 50000, xlim = c(15000,1000000))

Créer deux jeux de données (maisons / appartements)

table(etape5$type)

Appartement      Maison 
      42003       68734 
Maisons <- etape5 %>% filter(type == 'Maison')
Appartement <- etape5 %>% filter (type == 'Appartement')    

Fixer un seuil maximal des surfaces (histogramme)

par(mfrow=c(1,2))
hist(Maisons$surface, nclass = 500, xlim = c(0,600))
hist(Appartement$surface, nclass = 500, xlim = c(0,300))

Etape 6 : Sélection des bornes de prix et de surface

etape6 <- etape5 %>% 
       filter(between(prix, 15000, 5000000)) %>%
       filter(case_when(type=='Appartement' ~  
                          between(surface, 10, 200)) | 
       case_when(type=='Maison' ~ 
                   between(surface, 10, 350)))

Etape 7 : Création de la variable prix/m2

etape7 <- etape6 %>% mutate(prixm2 = prix/surface)

Etape 8 : Sélection des bornes de prix au m2

Fixer un seuil minimal des prix au m2 (percentile)

quantile(etape7$prixm2,c(0.01,0.99))
       1%       99% 
 387.3057 5825.8048 

Fixer un seuil maximal des prix au m2 (histogramme)

hist(etape7$prixm2, breaks = 1000, xlim = c(0,10000))

Filtrer les valeurs aberrantes

etape8 <- etape7 %>% 
          filter(between(prixm2,330,8000))

Etape 9 : Corrections diverses

Transformation de la date en année

etape8$date <- as.character(etape8$date)
etape8 <- etape8 %>% 
            mutate(annee = substr(etape8$date, 1, 4))

Arrondir les variables numériques

etape8$prix <- round(etape8$prix)
etape8$prixm2 <- round(etape8$prixm2)

Rétablir le ‘0’ du code communal et départmental

code <-as.character(etape8$INSEE_COM)
code[nchar(code)==4]<-paste0("0",code[nchar(code)==4])
etape8$INSEE_COM<-code
etape8$INSEE_DEP<-substr(etape8$INSEE_COM,1,2)

Etape 10 : Exportation

Mise en ordre des lignes et colonnes

DVFOK <- etape8 %>% select(id, date, annee = annee, 
                           INSEE_COM, INSEE_DEP,
                           type, prix, surface, prixm2,
                           latitude, longitude)%>%
                     arrange(date)

Ecrire le jeu de données final en csv

write.csv(DVFOK, 
          paste0(myrep, "dvf_clean.csv"),
          quote = FALSE)

Ecrire le jeu de données final en .RDS

saveRDS(DVFOK, paste0(myrep, "dvf_clean.RDS"))

Bilan

Fichier initial

dim(DVF)[1]
[1] 151749

Fichier nettoyé

dim(DVFOK)[1]
[1] 108695

% de lignes supprimées

100*(dim(DVF)[1]-dim(DVFOK)[1])/dim(DVF)[1]
[1] 28.37185

DVF3 : Cartographie

Retour sur sf

Nous revenons sur le package sf (spatial features) que nous avons déjà rencontré au moment de la création de cartes thématiques par IRIS ou communes à l’aide du package mapsf.

Ici le package sf va être utilisé pour cartographier non pas des zones mais des localisations ponctuelles. Il pourra être à nouveau couplé avec le logiciel de cartogaphie statique comme mapsf , afin par exemple de placer les localisations des ventes foncières au dessus du fonds de carte des IRIS ou communes.

Mais il pourra aussi servir de base à des cartographies dynamiques permettant de placer les points correspondant aux ventes sur des réseaux de rue et plus généralement sur des “tuiles” cartographiques permettant d’effectur des zoom. On utilisera à cet effet d’autres packages comme leaflet ou sa version simplifiée mapview.

Zone d’étude

On définit les paramètres de la zone d’étude en ajoutant le lien vers les dvf nettoyées

# Commune centre 
codectr<-  "35238"  # Code
namectr <- "Rennes"  # Nom
latctr <-  -1.67789
lonctr <-  48.11204
# Choix du rayon de collecte des dvf (en mètres)
rayon <- 50000
# Dossier de stockage
myrep <- "data/Rennes/"
# Fichier dvf nettoyé
mydvf <-"dvf_clean.RDS"

Choix de communes

On choisit plus précisément de travailler sur les communes de l’agglmomération

map_com <- readRDS(paste0(myrep,"EPCI_mapcom.RDS")) %>%
           select(INSEE_COM, NOM, geometry)
plot(map_com$geometry,
     col="lightyellow")

Sélection des dvf nettoyées

Nous reprenons le fichier nettoyé de localisation des ventes de maison et d’appartement établi au chapitre précédent et nous ne conservons que 9 variables ainsi que les ventes situées dans les communes de la zone choisie

dvf <- readRDS(paste0(myrep,mydvf)) %>%
        select(annee, INSEE_COM,
               type, prix, surface, prixm2,
               latitude, longitude) %>%
        filter(INSEE_COM %in% map_com$INSEE_COM)


kable(head(dvf,3))
annee INSEE_COM type prix surface prixm2 latitude longitude
2014 35238 Appartement 81000 66 1227 48.11526 -1.692243
2014 35238 Appartement 98500 33 2985 48.11689 -1.681832
2014 35238 Appartement 121000 62 1952 48.10974 -1.705489

Transformation des dvf en fichier sf

On ajoute une colonne geometry en utilisant la fonction st_as_sf()du package sf. On attribue à la projection le code EPSG = 4326 qui correspond à un fichier latitude longitude.

map_dvf <- st_as_sf(dvf, coords = c("longitude","latitude"))
st_crs(map_dvf)<- 4326

kable(head(map_dvf,3))
annee INSEE_COM type prix surface prixm2 geometry
2014 35238 Appartement 81000 66 1227 POINT (-1.692243 48.11526)
2014 35238 Appartement 98500 33 2985 POINT (-1.681832 48.11689)
2014 35238 Appartement 121000 62 1952 POINT (-1.705489 48.10974)

Test de superposition

On tente de superposer les deux cartes des ventes et des communes.

Commune puis dvf

plot(map_com$geometry, col="lightyellow")
plot(map_dvf$geometry, col="red", cex=0.2, pch=20,add=T)

dvf puis communes

plot(map_dvf$geometry, col="red", cex=0.2, pch=20)
plot(map_com$geometry, col="lightyellow", add=T)

  • Commentaire : les deux fonds de carte existent mais ils ne se superposent pas. Il y a donc un problème de projection !

Vérification de la projection

Nous savons que les coordonnées du fichier DVF sont en latitude longitude (EPSG=4326). Mais la projection de notre fonds des communes est en Lambert-93 (EPSG = 2154)

st_crs(map_com)
Coordinate Reference System:
  User input: RGF93 v1 / Lambert-93 
  wkt:
PROJCRS["RGF93 v1 / Lambert-93",
    BASEGEOGCRS["RGF93 v1",
        DATUM["Reseau Geodesique Francais 1993 v1",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4171]],
    CONVERSION["Lambert-93",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",46.5,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",3,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",49,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",44,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",700000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",6600000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["France - onshore and offshore, mainland and Corsica."],
        BBOX[41.15,-9.86,51.56,10.38]],
    ID["EPSG",2154]]

A priori il s’agit bien de la même de sorte que les coordonnées X,Y devraient bien se superposer sur le fonds IRIS

Ajustement des deux projections

On transforme donc le fichier des DVF en Lambert-93 pour obtenir une adéquation avec le fichier des communes.

map_dvf <- st_transform(map_dvf,2154)

Et cette fois-ci on réussit à superposer les deux fonds de carte

plot(map_com$geometry, col="lightyellow")
plot(map_dvf$geometry, col="red", cex=0.2, pch=20,add=T)

Cartographie ponctuelle avec mapsf

On peut utiliser les fonctions de mapsf vues précédemment. Par exemple repérer les ventes de maison et d’appartement :

mf_map(map_com, type="base",
       col="lightyellow")
mf_map(map_dvf, type = "typo",
       var= "type",add =T,
       border=NA,
       cex=0.3,
       leg_title = "Bien vendu")
mf_label(map_com,
         var="NOM",cex = 0.4,
         overlap = F,
         col="black",
         halo=T)
mf_layout(title = "Ventes de maisons et d'appartement 2013-2020",
          credits = "Sources : IGN, INSEE, DVF",
          scale = T)

Cartographie ponctuelle avec mapview

Le package mapview pemet de créer facilement des cartes interactives ù l’on pourra visualiser les ventes sur un fonds de carte de son choix (OSRM, ESRI, …). On peut l’utiliser par exemple pour voirune commune précise

library(mapview)
map_dvf_sel<-map_dvf %>% filter(INSEE_COM == "35196")
mapview(map_dvf_sel)

On peut améliorer sensiblement la visualisation à l’aide de paramètres présentés sur le site web du package https://r-spatial.github.io/mapview/

library(mapview)
# Choix de la zone d'étude
map_com_sel<-map_com %>% filter(INSEE_COM == "35196")
map_dvf_sel<-map_dvf %>% filter(INSEE_COM == "35196")

# Choix des tuiles
mapviewOptions(basemaps = c("OpenStreetMap" ,
                           "Esri.WorldImagery"))
# Première couche
map1 <- mapview(map_com_sel,
                col.regions="lightyellow",
                alpha.regions = 0.4)

# Deuxième couche
map2<-mapview(map_dvf_sel,
              zcol = "type",
              cex = 4,
              col.regions=c("blue","red"),
             alpha.regions = 0.4)
# Assemblage
map1+map2

Agrégation par communes

On peut évidemment agréger les informations sur les valeurs foncières par communes. Par exemple compter le nombre de maisons vendues et leur prix moyen au m2

dvf_by_com <- dvf %>% 
               group_by(INSEE_COM) %>%
               filter(type == "Maison") %>%
               summarise(nb = n(),
                         med_prixm2 = median(prixm2))

kable(head(dvf_by_com,3))
INSEE_COM nb med_prixm2
35001 400 2417.0
35022 66 1369.5
35024 745 2576.0

On peut désormais effectuer la jointure entre les données agrégées par adresse et le fichier sf de localisation des adresses :

map_com_dvf <- inner_join(map_com,dvf_by_com) %>% st_as_sf()

On peut désormais utiliser les méthodes de cartographie déjà vues avec mapsf :

mf_theme("agolalight")
mybreaks = c(1000, 1500, 2000, 2500, 3000, 3500,
             4000, 4500, 5000, 5500, 6000)
mypal=brewer.pal(n = 10,name = "Spectral")
mf_map(map_com_dvf, type = "choro", 
       var = "med_prixm2", 
       breaks=mybreaks, pal = mypal,
       border="gray", lwd=0.3,
       leg_title = "prix en €/m2", leg_val_rnd = 0,
       leg_pos = "left")

mf_map( map_com_dvf,  type = "prop",
        var = "nb", inches = 0.12, 
        col="black",border="white",
        leg_title = "nb. de ventes",
        leg_pos = "topleft")

mf_layout(title = "Les ventes de maison de 2013 à 2020",
        frame = TRUE, credits = "Sources : IGN et DVF", 
        arrow = F )

Cartographie par IRIS

Supposons maintenant que l’on souhaite cartographier les prix des maisons par IRIS à l’intérieur de la commune de Rennes. Comment faire dans la mesure où l’IRIS n’est pas mentionné dans le fichier DVF ? Réponse : en superposant les deux fonds de carte et en les intersectant.

Superposition

map_iris_ctr<-readRDS(paste0(myrep,"CTR_mapiris.RDS")) %>%
               select("NOM_COM","INSEE_COM","CODE_IRIS","NOM_IRIS","geometry")
map_dvf_ctr <- map_dvf %>% filter(INSEE_COM %in% map_iris_ctr$INSEE_COM)
plot(map_iris_ctr$geometry, col="lightyellow")
plot(map_dvf_ctr$geometry,col="red",pch=21,cex=0.1,add=T)

Intersection

On peut faire l’opération avec un SIG ou bien avec le package sf en utilisant sa fonction st_intersect() qui va ajouter tous les attributs du fichier IRIS au fichier DVF :

inter<-st_intersection(map_dvf_ctr,map_iris_ctr)
kable(head(inter,4))
annee INSEE_COM type prix surface prixm2 NOM_COM INSEE_COM.1 CODE_IRIS NOM_IRIS geometry
62 2014 35238 Appartement 80000 44 1818 Rennes 35238 352381204 Les Cloteaux POINT (350454.3 6786503)
130 2014 35238 Appartement 70000 55 1273 Rennes 35238 352381204 Les Cloteaux POINT (350613.6 6786410)
158 2014 35238 Appartement 115000 79 1456 Rennes 35238 352381204 Les Cloteaux POINT (350694.4 6786355)
212 2014 35238 Appartement 109000 73 1493 Rennes 35238 352381204 Les Cloteaux POINT (350662.4 6786511)

Agrégation

Il ne reste plus qu’à refaire l’agrégation en suivant la même procédure que pour les communes mais en retirant les iris où il y a eu moins de 5 ventes.

map_iris_dvf <- inter %>% 
               st_drop_geometry() %>%
               group_by(CODE_IRIS, NOM_IRIS) %>%
               filter(type == "Maison") %>%
               summarise(nb = n(),
                         med_prixm2 = median(prixm2)) %>%
              filter(nb >=5) %>% 
              right_join(map_iris_ctr) %>% 
              st_as_sf()
               
             
kable(head(map_iris_dvf,3))
CODE_IRIS NOM_IRIS nb med_prixm2 NOM_COM INSEE_COM geometry
352380101 Cathédrale 5 2680 Rennes 35238 MULTIPOLYGON (((351427 6789…
352380102 Hoche 15 4671 Rennes 35238 MULTIPOLYGON (((352032.6 67…
352380109 Vieux Saint-Étienne 17 5263 Rennes 35238 MULTIPOLYGON (((351501.3 67…

Cartographie

Il ne reste plus qu’à faire la carte en reprenant le programme utilisé pour les communes :

mf_theme("agolalight")
mybreaks = c(1000, 1500, 2000, 2500, 3000, 3500,
             4000, 4500, 5000, 5500, 6000)
mypal=brewer.pal(n = 10,name = "Spectral")
mf_map(map_iris_dvf, type = "choro", 
       var = "med_prixm2", 
       breaks=mybreaks, pal = mypal,
       border="gray", lwd=0.3,
       leg_title = "prix en €/m2", leg_val_rnd = 0,
       leg_pos = "left")

mf_map(map_iris_dvf,  type = "prop",
        var = "nb", inches = 0.08, 
        col="black",border="white",
        leg_title = "nb. de ventes",
        leg_pos = "topleft")

mf_layout(title = "Les ventes de maison de 2013 à 2020",
        frame = TRUE, credits = "Sources : IGN et DVF", 
        arrow = F )

MOD1 : Analyse sociologique / Chi2

On se propose dans un premier temps de modéliser des caractéristiques non-spatiales des individus à partir des données du RP 2019. On va considérer ici uniquement des attributs catégoriels ce qui nous conduira à formuler des hypothèses aboutissant à la création d’un tableau de contingence et la réalisation d’un test du chi-2.

Définir le sujet

Soit le sujet : Qui sont les nouveaux propriétaires ? Quel est leur niveau de qualification ?

Définir le statut d’occupation

Propriétaire ? Locataire ? Occupant à titre gratuit ?

Définir la notion de “qualification” ?

Le diplôme le plus élevé ? le nombre d’années d’étude ?

Définir la date

Année 2019 uniquement ? Résultats du RP 2019 (2017-2021) ?

Définir la population cible

Personnes installées récemment ? depuis quand ?

Définir la zone d’étude

Commune centre ? Agglomération ? Aire urbaine ? …

Formuler des hypothèses

Justes ou fausses, les hypothèses permettent de cadrer l’analyse. Elles peuvent faire intervenir plusieurs variables

Diplôme et propriété

Les nouveaux propriétaires sont plus souvent diplômés.

Âge et propriété

Les propriétaires sont plus âgés que les locataires.

Propriété et territoire

Les propriétaires sont plus nombreux en zone rurale

Propriété, âge, diplôme, territoire

Les jeunes ménages sont locataires en ville avant de devenir propriétaires dans les zones périurbaines ou rurales. Ils sont plus ou moins proches du centre selon leur niveau de qualification qui détermine à long terme leur revenu …

Charger les données

On charge ici les données ménages de l’aire urbaine.

tab_ind<-readRDS(paste0(myrep,"AA_men.RDS"))
head(tab_ind[,1:5],2)
FALSE    COMMUNE   ARM      IRIS ACHL AEMM
FALSE 1:   22036 ZZZZZ ZZZZZZZZZ  A12 1973
FALSE 2:   22036 ZZZZZ ZZZZZZZZZ  B11 1992

Préparer l’analyse

  • Soit la relation entre statut d’occupation (Y) et diplôme le plus élevé du chef de ménage (X). Il s’agit de deux variables catégorielles (= qualitatives) que l’on va typiquement mettre en relation à l’aide d’un tableau de contingence et d’un test du chi-2. L’analyse statistique est simple sous R mais il faut tenir compte de trois difficultés

  • Le choix de la population de référence est important. Ici on va sélectionner les ménages installés depuis moins de 3 ans

  • la sélection ou le regroupement des diplômes est également important car cela va influer sur les résultats du test.

  • la pondération des individus doit également être prise en compte puisque le recensement est basé sur un sondage dans les zones urbaines

Sélection des individus et des variables

tab_sel<- tab_ind %>% 
  filter(as.numeric(ANEM) < 3) %>%
  select(INSEE_COM, ANEM, DIPLM,STOCD, IPONDL) 


knitr::kable(head(tab_sel,4))
INSEE_COM ANEM DIPLM STOCD IPONDL
22036 2 15 10 1.034483
22036 2 13 23 1.034483
22036 2 13 21 1.034483
22036 1 13 10 1.034483

Recodage des modalités

On cherche le code des modalités dans le fichier des métadonnées

meta<-readRDS(paste0(myrep,"men_meta.RDS"))
metasel <- meta %>% 
  filter(COD_VAR %in% c("STOCD", "DIPLM"))
kable(metasel[,c(1,3,4)])
COD_VAR COD_MOD LIB_MOD
DIPLM 01 Pas de scolarité ou arrêt avant la fin du primaire
DIPLM 02 Aucun diplôme et scolarité interrompue à la fin du primaire ou avant la fin du collège
DIPLM 03 Aucun diplôme et scolarité jusqu’à la fin du collège ou au-delà
DIPLM 11 CEP (certificat d’études primaires)
DIPLM 12 BEPC, brevet élémentaire, brevet des collèges, DNB
DIPLM 13 CAP, BEP ou diplôme de niveau équivalent
DIPLM 14 Baccalauréat général ou technologique, brevet supérieur, capacité en droit, DAEU, ESEU
DIPLM 15 Baccalauréat professionnel, brevet professionnel, de technicien ou d’enseignement, diplôme équivalent
DIPLM 16 BTS, DUT, Deug, Deust, diplôme de la santé ou du social de niveau bac+2, diplôme équivalent
DIPLM 17 Licence, licence pro, maîtrise, diplôme équivalent de niveau bac+3 ou bac+4
DIPLM 18 Master, DEA, DESS, diplôme grande école niveau bac+5, doctorat de santé
DIPLM 19 Doctorat de recherche (hors santé)
DIPLM YY Hors résidence principale
STOCD 00 Logement ordinaire inoccupé
STOCD 10 Propriétaire
STOCD 21 Locataire ou sous-locataire d’un logement loué vide non HLM
STOCD 22 Locataire ou sous-locataire d’un logement loué vide HLM
STOCD 23 Locataire ou sous-locataire d’un logement loué meublé ou d’une chambre d’hôtel
STOCD 30 Logé gratuitement

On recode les modalités des deux variables en regroupant certaines.

tab_sel$STOCD<-as.factor(tab_sel$STOCD)
levels(tab_sel$STOCD)<-c("Proprétaire","Locataire",
                         "Locataire","Locataire",NA)
tab_sel$DIPLM<-as.factor(tab_sel$DIPLM)
levels(tab_sel$DIPLM) <- c("Aucun","Aucun","Aucun",
                           "BEP","BEP","BEP", "BAC","BAC",
                         "BAC+123","BAC+123","> BAC+3","> BAC+3",NA)

                        
knitr::kable(head(tab_sel,3))
INSEE_COM ANEM DIPLM STOCD IPONDL
22036 2 BAC Proprétaire 1.034483
22036 2 BEP Locataire 1.034483
22036 2 BEP Locataire 1.034483

Tableau de contingence (FAUX)

Le plus simple semble être l’instruction table()

tab_cont<-table(tab_sel$STOCD,tab_sel$DIPLM)

knitr::kable(addmargins(tab_cont))
Aucun BEP BAC BAC+123 > BAC+3 Sum
Proprétaire 637 4132 3500 5738 3305 17312
Locataire 2944 9422 13158 12589 6237 44350
Sum 3581 13554 16658 18327 9542 61662

Tableau de contingence (JUSTE)

On pondère avec wtd.table() du package questionr.

programme

library(questionr)
tab_cont_wtd<-wtd.table(tab_sel$STOCD,tab_sel$DIPLM,
                        weights = tab_sel$IPONDL)

knitr::kable(round(addmargins(tab_cont_wtd),0))
Aucun BEP BAC BAC+123 > BAC+3 Sum
Proprétaire 751 4779 4167 7479 5125 22301
Locataire 4133 12437 21622 20839 11420 70450
Sum 4884 17216 25789 28317 16545 92751

Comparaison des résultats

  • Tableau non pondéré … légèrement faux !
Aucun BEP BAC BAC+123 > BAC+3 All
Proprétaire 17.8 30.5 21 31.3 34.6 28.1
Locataire 82.2 69.5 79 68.7 65.4 71.9
Total 100.0 100.0 100 100.0 100.0 100.0
  • Tableau pondéré … juste !
Aucun BEP BAC BAC+123 > BAC+3 All
Proprétaire 15.4 27.8 16.2 26.4 31 24
Locataire 84.6 72.2 83.8 73.6 69 76
Total 100.0 100.0 100.0 100.0 100 100

Visualisation du tableau de contingence

On choisit l’orientation du tableau et on l’affiche avec plot()

mytable<-wtd.table(tab_sel$DIPLM,tab_sel$STOCD,weights = tab_sel$IPONDL)
plot(mytable)

Tant qu’à faire, on améliore la figure avec des paramètres supplémentaires :

plot(mytable, 
     main = "Propriété et diplôme dans l'aire urbaine", 
     sub = "Source : INSEE - RP 2019 - Ménages installés depuis - de 3 ans",  
     col=c("lightyellow","lightgreen"))

Test du Chi-deux

Ce test se réalise facilement sur le tableau de contingence avec l’instruction chisq.test() :

mytest<-chisq.test(mytable)
mytest

    Pearson's Chi-squared test

data:  mytable
X-squared = 1731.7, df = 4, p-value < 2.2e-16
  • Interprétation : La relation est très significative (p<0.001)

Prolongements

On peut introduire toute une série de variantes à partir de ce premier programme. Par exemple, si on retire la commune-centre, les résultats sont assez différents :

tab_sel2<- tab_sel %>% 
  filter(INSEE_COM != codectr)

mytable<-wtd.table(tab_sel2$DIPLM,tab_sel2$STOCD,weights = tab_sel2$IPONDL)

plot(mytable, 
     main = "Propriété et diplôme hors commune-centre", 
     sub = "Source : INSEE - RP 2019 - Ménages installés depuis - de 3 ans",  
     col=c("lightyellow","lightgreen"))

chisq.test(mytable)

    Pearson's Chi-squared test

data:  mytable
X-squared = 733.11, df = 4, p-value < 2.2e-16

A vous de deviner pourquoi :-)


MOD2 : Analyse spatiale / Gradient

On se propose d’analyser la variation d’un phénomène quantitatif (le prix de vente des maisons en € par m2) en fonction de la distance à un point remarquable (le centre de l’agglomération). On fait l’hypothèse que ce prix correspond à un gradient de décroissance, conformément aux théories de la rente foncière. Pour modéliser ce gradient, nous ferons appels à un modèle de régression linéaire ou non-linéaire.

Formuler l’hypothèse

Soit l’ensemble des ventes effectuées dans un rayon de 50 km autour d’une agglomération. Existe-t-il une relation entre (Y) le prix de vente mesuré en € par m2 de surface habitable et (X) la distance au point central de l’agglomération ?

N.B. Nous retirons la commune centre de l’analyse car on suppose que son comportement interne obéit à une logique différente

Définir l’espace d’étude

Nous allons retenir la distance à vol d’oiseau entre chaque point ou a eu lieu une vente et le centre de l’agglomération que l’on fixera en fonction de notre connaissance du terrain.

# Commune centre 
codectr<-  "35238"  # Code
namectr <- "Rennes"  # Nom

# Choix du rayon de collecte des dvf (en mètres)
rayon <- 50000

# Dossier de stockage
myrep <- "data/Rennes/"

# Fichier dvf nettoyé
mydvf <-"dvf_clean.RDS"

Calculer les distances au centre

On doit ajouter à notre fichier dvf la distance au point central retenu.

## Fichier des dvf
dvf<-readRDS(paste0(myrep,mydvf))

# Aire urbaine
map_com<-readRDS(paste0(myrep,"AA_mapcom.RDS"))

# Selection des dvf
dvf<-dvf %>% filter(INSEE_COM %in% map_com$INSEE_COM) %>%
             filter(INSEE_COM != codectr)

# Transformation cartographique des dvf
map_dvf <- st_as_sf(dvf,coords =c("longitude","latitude")) 
st_crs(map_dvf)<-4326
map_dvf<-st_transform(map_dvf,2154)

## Commune centre => centroïde
map_ctr <-readRDS(paste0(myrep,"CTR_mapcom.RDS"))
map_ctr <-st_centroid(map_ctr)

## Calcul de la distance
map_dvf$dist<-as.numeric(st_distance(map_dvf,map_ctr))/1000

## Visualisation rapide
plot(map_dvf[,"dist"])

Modèle linéaire

On teste l’existence d’un modèle linéaire en fixant les paramètres suivants :

Paramètres

Y <-map_dvf$prixm2
X <-map_dvf$dist
nomY <- "Prix au m2"
nomX <- "Distance en km"
model<- "Modèle Linéaire : Y = aX+b"

Ajustement

mod1 <- lm(Y~X)
summary(mod1)

Call:
lm(formula = Y ~ X)

Residuals:
    Min      1Q  Median      3Q     Max 
-2309.1  -368.5   -47.0   314.5  5856.8 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2837.7774     5.9202   479.3   <2e-16 ***
X            -45.1457     0.3075  -146.8   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 633.1 on 50618 degrees of freedom
Multiple R-squared:  0.2986,    Adjusted R-squared:  0.2986 
F-statistic: 2.155e+04 on 1 and 50618 DF,  p-value: < 2.2e-16

Visualisation

plot(X,Y,
     xlab = nomX,
     ylab = nomY,
     main = model,
     pch=20,
     col="red",
     cex=0.1,
     sub = " Source : DVF 2013-2020"
     )

abline(mod1, col="blue", lwd=2)

Modèle exponentiel

On teste l’existence d’un modèle exponetielle en fixant les paramètres suivants :

Paramètres

Y <-log(map_dvf$prixm2)
X <-map_dvf$dist
nomY <- "log(Prix au m2)"
nomX <- "Distance en km"
model<- "Modèle Exponentiel : log(Y) = aX+b"

Ajustement

mod2 <- lm(Y~X)
summary(mod2)

Call:
lm(formula = Y ~ X)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.95482 -0.14476  0.02937  0.19630  1.68254 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.9791861  0.0031466  2535.8   <2e-16 ***
X           -0.0245287  0.0001635  -150.1   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3365 on 50618 degrees of freedom
Multiple R-squared:  0.3079,    Adjusted R-squared:  0.3079 
F-statistic: 2.252e+04 on 1 and 50618 DF,  p-value: < 2.2e-16

Visualisation

plot(X,Y,
     xlab = nomX,
     ylab = nomY,
     main = model,
     pch=20,
     col="red",
     cex=0.1,
     sub = " Source : DVF 2013-2020"
     )

abline(mod2, col="blue", lwd=2)

Modèle logarithmique

On teste l’existence d’un modèle logarithmique en fixant les paramètres suivants :

Paramètres

Y <-map_dvf$prixm2
X <-log(map_dvf$dist)
nomY <- "Prix au m2"
nomX <- "log(Distance en km)"
model<- "Modèle Logarithmique : Y = a.log(X)+b"

Ajustement

mod3 <- lm(Y~X)
summary(mod3)

Call:
lm(formula = Y ~ X)

Residuals:
    Min      1Q  Median      3Q     Max 
-2820.5  -364.0   -31.8   320.7  5900.1 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 3801.054     12.449   305.3   <2e-16 ***
X           -649.871      4.559  -142.5   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 638.5 on 50618 degrees of freedom
Multiple R-squared:  0.2864,    Adjusted R-squared:  0.2864 
F-statistic: 2.032e+04 on 1 and 50618 DF,  p-value: < 2.2e-16

Visualisation

plot(X,Y,
     xlab = nomX,
     ylab = nomY,
     main = model,
     pch=20,
     col="red",
     cex=0.1,
     sub = " Source : DVF 2013-2020"
     )

abline(mod3, col="blue", lwd=2)

Modèle puissance

On teste l’existence d’un modèle puissance en fixant les paramètres suivants :

Paramètres

Y <-log(map_dvf$prixm2)
X <-log(map_dvf$dist)
nomY <- "log(Prix au m2)"
nomX <- "log(Distance en km)"
model<- "Modèle Logarithmique : Y = a.log(X)+b"

Ajustement

mod4 <- lm(Y~X)
summary(mod4)

Call:
lm(formula = Y ~ X)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.13041 -0.14325  0.04097  0.20474  1.55332 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.467618   0.006720  1260.1   <2e-16 ***
X           -0.339949   0.002461  -138.1   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3447 on 50618 degrees of freedom
Multiple R-squared:  0.2738,    Adjusted R-squared:  0.2738 
F-statistic: 1.908e+04 on 1 and 50618 DF,  p-value: < 2.2e-16

Visualisation

plot(X,Y,
     xlab = nomX,
     ylab = nomY,
     main = model,
     pch=20,
     col="red",
     cex=0.1,
     sub = " Source : DVF 2013-2020"
     )

abline(mod4, col="blue", lwd=2)

Choix du meilleur modèle

On décide ici de retenir le modèle 2 (exponentiel) et de stocker les valeurs estimées et résiduelles des prix au m2.

mod<-mod3
map_dvf$prixm2_estim <- mod$fitted.values
map_dvf$prixm2_resid <- mod$residuals

Agrégation par commune

On agrège les valeurs estimées et résiduelles par commune

mod_com <- st_drop_geometry(map_dvf) %>%
              group_by(INSEE_COM) %>%
              summarise(nb=n(),
                        prixm2 = mean(prixm2),
                        prixm2_estim = mean(prixm2_estim),
                        prixm2_resid = mean(prixm2_resid)) %>%
              arrange(prixm2_resid)

On en déduit les communes où le prix est plus bas que ce que laisserait prévoir la distance au centre :

kable(head(mod_com), digits=0)
INSEE_COM nb prixm2 prixm2_estim prixm2_resid
35281 1878 2435 2999 -564
35261 18 1085 1571 -486
35242 37 1045 1522 -476
35244 31 1093 1535 -443
35290 75 1122 1554 -432
35164 80 1164 1577 -413

Et celles où le prix est plus élevé que ce que laisserait prévoir la distance au centre

kable(tail(mod_com),digits=0)
INSEE_COM nb prixm2 prixm2_estim prixm2_resid
35001 655 2401 2199 201
35069 1059 2225 2016 209
35177 486 2383 2162 222
35152 766 2200 1956 243
35047 2012 2533 2270 263
35051 1560 3134 2720 414

Cartographie des résultats

En se limitant à l’aire urbaine, on va visualiser les résultats du modèle. On commence par effectuer la jointure entre les résultats du modèle de régression et le fonds de carte des communes :

map_com<-map_com %>%
         left_join(mod_com)

On peut alors cartographier les prix observés et théoriques :

mybreaks <-c(1000,1500, 1750, 2000, 2250, 2500, 2750,  3000, 6000)
mypal <- brewer.pal(n = 8, name="YlOrRd")

par(mfrow = c(1,2))
mf_map(map_com, 
       type = "choro",
       var="prixm2_estim",
       breaks=mybreaks,
       pal=mypal,
       leg_val_rnd = 0,
       col_na = "gray80",
       leg_pos = "topleft")
mf_layout("Prix théoriques", 
          arrow=F, frame=T,
          credits = "Source : DVF, 2013-2020")


mf_map(map_com, 
       type = "choro",
       var="prixm2",
       breaks=mybreaks,
       pal=mypal,
       leg_val_rnd = 0,
        col_na = "gray80",
       leg_pos = "topleft")
mf_layout("Prix observés",
          arrow=F, frame=T,
          credits = "Source : DVF, 2013-2020")

Et en déduire la carte des résidus :

mybreaks <-c(-700, -500, -300,-100,100,300,500,700)
mypal <- brewer.pal(n = 7, name="RdBu")
par(mfrow = c(1,1))
mf_map(map_com, 
       type = "choro",
       var="prixm2_resid",
       breaks=mybreaks,
       pal=mypal,
       col_na = "gray80",
       leg_val_rnd = 0,
       leg_pos = "topleft")
mf_layout("Prix résiduels", 
          arrow=F, frame=T,
          credits = "Source : DVF, 2013-2020")